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ERROR ANALYSIS OF TAU-LEAP SIMULATION METHODS 

By David F. Anderson * Arnab Ganguly * and Thomas G. Kurtz * 

University of Wisconsin - Madison* 

We perform an error analysis for numerical approximation methods of 
continuous time Markov chain models commonly found in the chemistry and 
biochemistry literature. The motivation for the analysis is to be able to com- 
pare the accuracy of different approximation methods and, specifically, Euler 
tau-leaping and midpoint tau-leaping. We perform our analysis under a scal- 
ing in which the size of the time discretization is inversely proportional to 
some (bounded) power of the norm of the state of the system. We argue that 
this is a more appropriate scaling than that found in previous error analyses in 
which the size of the time discretization goes to zero independent of the rest 
of the model. Under the present scaling we show that midpoint tau-leaping 
achieves a higher order of accuracy, in both a weak and a strong sense, than 
Euler tau-leaping; a result that is in contrast to previous analyses. We present 
examples that demonstrate our findings. 

1. Introduction. This paper provides an error analysis for numerical approximation methods for con- 
tinuous time Markov chain models that are becoming increasingly common in the chemistry and biochem- 
istry literature. Our goals of the paper are two-fold. First, we want to demonstrate the importance of consid- 
ering appropriate scalings in which to carry out error analyses for the methods of interest. Second, we wish 
to provide such an error analysis in order to compare the accuracy of two different approximation methods. 
We perform our analysis on the Euler tau-leaping method first presented in [11] and a midpoint tau-leaping 
method developed below, which is only a slight variant of one presented in [11]. The midpoint tau-leaping 
method will be demonstrated to be more accurate than Euler tau-leaping in both a strong and a weak sense, a 
result that is in contrast to previous error analyses. We will discuss why previous error analyses made differ- 
ing predictions than does ours and argue that the scaling provided here, or variants thereof, is a more natural 
and appropriate choice for error analyses of such methods. We also provide examples that demonstrate our 
findings. 

1.1. The basic model. The motivation for the class of mathematical models under consideration comes 
from chemistry and biochemistry, and more generally from population processes (though we choose the lan- 
guage of chemistry throughout the paper). We assume the existence of a chemical reaction system consisting 
of (i) d chemical species {5i, 52, . . . , Sd} and (ii) a finite set of possible reactions, which we index by k. 
Each reaction requires some number of the species as inputs and provides some number of the species as 
outputs. For example, the reaction Si — )• 2^2 would require one molecule of Si for the input and provide 
two molecules of S2 for the output. If reaction k occurs at time t, then the state of the system X{t) E Z>q 
is updated via addition of the reaction vector £ ^'^> which represents the net change in the abundances 
of the underlying species: 

X{t) = X{t-) + Uk- 

Returning briefly to the example Si — )■ 252, the associated reaction vector for this reaction would be 
[—1, 2,0,..., 0]^. Finally, we denote by vf. the vector in Z>q representing the source of the fcth reaction. 
Returning again to the example Si — 252, the source vector for this reaction is = [1, 0, . . . , 0]^. 
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We assume that the waiting times for the k reactions are exponentially distributed with intensity functions 
Afc : M>Q — M>o. We extend each Afc to all of by setting it to zero outside M>q. This model is a 
continuous time Markov chain in Z>q with generator 

(1.1) {Af){x) = ^ Xk{x)if{x + uk) - fix)), 

k 

where / : Z'' — M is arbitrary. Kolmogorov's forward equation for this model, termed the "chemical master 
equation" in the chemistry and biology literature, is 

^P{x,t\-K) = '^P{x - i^k,t\TT)Xkix - Uk) - y^P(x,t|7r)Afc(x), 

k k 

where for x G Z>q P{x, t\Tr) represents the probability that X{t) = x, conditioned upon the initial distri- 
bution vr. One representation for path-wise solutions to this model uses a random time change of Poisson 
processes: 

(1.2) X{t) = X{0) + E >^k{X{s))ds^ Uk, 

where the Yk are independent, unit-rate Poisson processes (see, for example, [16]). Note that X{t) = X{t) — 
X^fc So ^k{X{s))ds Vk is a martingale with quadratic covariation matrix [X]t = Yl,k ^/q XkiX{s))ds^ ^k^k- 

A common choice of intensity function for chemical reaction systems, and the one we adopt throughout, 
is mass action kinetics. Under mass action kinetics, the intensity function for the fcth reaction is 

X \ _ Ck yr Xt\ 5l£.f "TT T 

where Ck is a positive constant and Ck is defined by the above equation. Mass action kinetics arises by 
thinking of c^At as the approximate probability that a particular set of the molecules needed in the A;th 
reaction will react over a time-period of size At, and then counting the number of ways such a reaction 
could happen. Implicit in the assumption of mass action kinetics is that the vessel under consideration is 
"well-stirred." For ease of notation we will henceforth drop the indicator functions from our representation 
of mass action kinetics. More general rates will be discussed in the remark at the top of page six. 

1 .2. Numerical methods. There are a number of numerical methods that produce statistically exact sam- 
ple paths for the model described above. These include the stochastic simulation algorithm, better known as 
Gillespie's algorithm ([9, 10]), the first reaction method ([9]), and the next reaction method ([1, 8]). All such 
algorithms perform the same two basic steps multiple times until a sample path is produced over a desired 
time interval: first, conditioned on the current state of the system the amount of time that passes until the 
next reaction takes place. A, is computed and second the specific reaction that has taken place is found. If, 
however, Ylik Xk{X{t)) S> then A ^ Xk{X{t)))~^ <^ 1 and the time needed to produce a single 
exact sample path over a time interval can be prohibitive. 

The approximate algorithm "tau-leaping" was developed by Dan Gillespie in [1 1] in an effort to overcome 
the problem that A may be prohibitively small. The basic idea of tau-leaping is to hold the intensity functions 
fixed over the time interval [tn,tn + h] at the values Xk{X{tn)), where X{tn) is the current state of the 
system, and, under this assumption, compute the number of times each reaction takes place over this period. 
As the waiting times for the reactions are exponentially distributed this leads to the following algorithm. 



(1.3) 



Xk{x) = Ck 
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Algorithm 1 (Euler tau-leaping). Set Z(0) = X{0), to = 0, n = and repeat the following until 

tn+l > T. 

1. Set Z{tn+i) = Z{tn) + J2k'Pk,n{^k{Z{tn))h)uk, Set tn+l = tn + h, and set n = n + 1, where 
Vk,n{x) are independent Poisson random variables with parameters x. 

Several improvements and modifications have been made to the basic algorithm described above over the 
years. However, they are mainly concerned with how to choose the step-size adaptively [4, 12] and/or how 
to ensure that population values do not go negative during the course of a simulation [2, 3, 5], and are not 
explicitly relevant to the current discussion. 

Similar to (1.2), a path-wise representation of Euler tau-leaping can be given through a random time 
change of Poisson processes: 



(1.4) Z{t) = X{0) + Y.Yk (^J^ Xk{Z o i]{s))ds^ 



T^k, 



where r]{s) = tn if tn < s < tn+i and the Yk are as before. Noting that /q"^^ Xk{Z o ri{s))ds = 
Y17=Q ^k{Z {ti)){ti+i — ti) explains our choice to call this method "Euler tau-leaping." Defining the op- 
erator 

(1.5) {B.f){x) = ^k{z){f{x + uk) - fix)), 

k 

we see that for t > 

(1.6) = Ef{Z o r]{t)) +E f {Bzo^it)f){Z{s))ds, 

Jr,{t) 

SO long as the expectations exist. Further, we note that Z{t) = Z{t) — lo ^k{Z o r]{s))dsvk is a 
martingale with quadratic covariation matrix [Z]t = J2k ^k ( Jq ^k{Z o r]{s))ds] Vk^k- 



It is natural to believe that a midpoint type method would be more accurate than an Euler type method in 
many situations. We therefore define the function 



p{z) = z + l^hy^Xk{z)uk, 



2 

k 

which computes an approximate midpoint for the system assuming the state of the system is z and the 
time-step is h. 

Algorithm 2 (Midpoint tau-leaping). Set Z{0) = X{0), to = 0,n = and repeat the following until 

tn+l > T. 

1. Set Z{tn+i) = Z{tn) + J2k 'Pk,n{^k ° P° Z{tn)h)i'k, sct tn+l = tn + h, and sct 71 = n + 1, where 
Vk,n{x) are independent Poisson random variables with parameters x. 

Similar to (1.2) and (1.4), Z{t) can be represented via a random time change of Poisson processes: 



Z{t)=X{^) + Y,Yk(^j'^ 



Xk o p o Z o r]{s)ds I Vk 
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where 7]{-) is as before. For Bz defined via (1.5) and any < t and any function / 

(1.7) = Ef{Z o rj{t)) + E f {BpoZo^^t)f){Z{s))ds. 

hit) 

Finally, Z{t) = Z{t) — J2k j^^k ° P ° Z, o rj{s)ds Uk is a martingale with quadratic covariation matrix 
[Z]t = Ylk (yj^ po Z o r]{s)ds^ ^k^k- The main goal of this paper is to show that the midpoint 
tau-leaping algorithm is indeed more accurate than the Euler tau-leaping method under an appropriate, and 
natural, scaling described in Section 2. 

Remark. Historically the time discretization parameter for tau-leaping has been r, thus giving the 
method its name. We choose to break from this tradition and denote our time-step by h so as not to confuse 
r with a stopping time. 

1.3. Previous error analyses. Under the scaling /i — Rathinam et al. performed a consistency check 
of Euler tau-leaping and found that the local truncation error was 0{h?) for all moments [18]. They also 
showed that under this same scaling Euler tau-leaping is first order accurate in a weak sense in the special 
case that the intensity functions Afc are linear [18]. Li extended these results by showing that as /i — 
Euler tau-leaping has a strong error (in the norm) of order 1/2 and a weak error of order one [17], which 
agree with classical results pertaining to numerical analysis of SDEs driven by Brownian motions (see, for 
example, [13]). 

Under the scaling /i — )• it is readily seen that midpoint tau-leaping is no more accurate than Euler 
tau-leaping. This follows since midpoint tau-leaping consists of making an 0{h?') correction to the intensity 
functions used in Euler tau-leaping. As /i — 0, this correction becomes negligible as Poisson processes 
"ignore" 0{h?) corrections, and the accuracy of the two methods will be the same. 

We simply note that while the analyses performed in [18] and [17] and the argument made in the previous 
paragraph are technically correct, performing an analysis as /i — 0, independent of the rest of the model, 
is at odds with the useful regime of tau-leaping. That is, tau-leaping would only be used in a regime where 
/i 3> A, where A is the expected amount of time between reactions, for otherwise an exact method would 
be performed. Therefore, we should require that 

(1.8) /i» (J]Afc(Z(t)))-i or /i^Afc(Z(t))»l, 

k k 

where Z{t) is the state of the system. In Section 2 we will present a natural scaling for the models under 
consideration that does satisfy (1.8) and under which we will perform our analysis. 

1.4. Paper outline. The remainder of the paper is organized as follows. In Section 2 we give some nat- 
ural assumptions on the models considered in this paper and introduce the scaling under which we perform 
our analysis. In Section 3 we perform a strong error analysis for both the Euler and midpoint tau-leaping 
methods and show that midpoint tau-leaping is the more accurate of the two under our scaling. In Section 4 
we perform a weak error analysis of the different methods and again conclude that the midpoint method is 
more accurate. In Section 5 we present numerical examples demonstrating our results. 

2. Assumptions on the model. 

2.1. Scalings of the model and the algorithms. As discussed in the introduction, tau-leaping methods 
will only be of use if the time-discretization parameter h satisfies h Y^j^ Xk{Z{t)) ^ 1 while Q2k ^kiZ{t)))~ 
1, where Z{t) is the state of the system at time t. There are a number of ways for the second condition to 
hold and a modeling choice must be made. We make the following natural assumptions: 
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(i) The initial abundance of each species scales with V for some V ^ 1. 

(a) Eachrate constant satisfies = ©(^^"^fe'^), where 1 = [1, 1, ... , 1]-^. In particular, = dk/V^~'^i' 
for some dk > 0. 

We will denote by the normalized process defined as the vector of abundances divided by V, and will 
denote by the intensity function defined to be mass action kinetics with rate constants . This scaling is 
the so called "classical scaling" and arises naturally by thinking of V as the volume of the vessel in which the 
reactions are taking place multiplied by Avogadro's number ([14]). In this case gives the concentration 
of each species in moles per unit volume. To understand the scaling for the rate constants, consider the 
case of a reaction requiring as input two constituent molecules: and 5*2. Perhaps Si + S2 ^ S3. It is 
reasonable to assume that the probability that a particular pair of Si and S2 molecules meet, and react, in a 
small time interval is inversely proportional to the volume of the vessel. This same type of logic holds for 
the cases in which more than two molecules are needed for a reaction to take place (i.e. the probability that 
three particular molecules meet and react is inversely proportional to the volume squared). For the case that 
only one molecule is needed for a reaction to take place, it is reasonable to assume that the probability of 
such a reaction taking place in the next small interval of time for a particular molecule should not scale with 
the volume. See also [19], Chapter 6. 

Models that satisfy assumptions {i) and {ii) above have an important property that we will detail here 
and make use of later. Let x{t) denote the solution to the deterministic initial value problem 



(2.1) x{t) = F{x{t)) = Y,dkx{tY^Uk, a;(0) = xoG 



■pd 



where is defined in assumption {ii) above, and where for any two vectors u" = ■ ■ ■ u'j'- and we adopt 
the convention that 0*^ = 1. That is, x{t) is the solution to the corresponding deterministically modeled 
chemical reaction system with mass action kinetics. It was shown in [14, 15] that for any e > and any 

T > 0, if X^(0) = a;(0) = xq, then 

(2.2) lim P{ sup \X^{t) - x{t)\ > e} ^ 0. 

V^oo fe[o,T] 

Denoting Ajt as deterministic mass action kinetics with rate constant d^, it is an exercise to check that for 
any reaction, i.e. zeroth order, first order, second order, etc., and any x G R>q 

Xl{Vx) = V\k{x) + Ql{x), 

where is uniformly bounded in V and is nonzero only if the reaction requires more than one molecule of 
a particular species as an input. For example, for the second order reaction 5i + 5*2 — )■ 5*3 we have 

Afc (^x) = ^ {Vxi) {Vx2) = VdkXiX2 = V\k{x), 
whereas for the second order reaction — ^ 5*3 we have 

>^l{yx) = '^Vxi {Vxi - 1) = Vdkxi - dkxi = V\k{x) + Ck{x). 



with C,^ {x) = —dkXi. The term will have a true V dependence if three or more molecules of a particular 

^>0 



species are required as input. We now make the definition (x) = y (Vx) , and note that for all x G 



(2.3) Al{x) = Xk{x) + ^Ckix), 

and A^(x) = if X ^ M>o- Manipulating the definition of shows that for all x eW^ 

(2.4) XliVx) = VAl{x). 
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Remark. The assumption of mass action kinetics is not critical to the analysis carried out in this 
paper. Instead, what is critical to this particular analysis is that our kinetics satisfies the scaling (2.4) for 
satisfying (2.3) with sufficiently smooth. 

We now choose a discretization parameter for the approximate methods that is dependent upon the as- 
sumptions of the model set out above. We let 



With this choice of time-step, we let Z and Z denote the processes generated by Euler and midpoint 
tau-leaping, respectively, normalized by V. We can now state more clearly what the analysis of this paper 
will entail. We will consider the case of y » by letting V ^ oo and consider the relationship of the 
normalized approximate processes and to the original process , normalized similarly. Note that 
all three processes converge to the solution of (2.1). We will perform both weak and strong error analyses. 
In the strong error analysis, we will consider convergence as opposed to the more standard (at least for 
systems driven by Brownian motions) convergence. The reason for this is simple: the Ito isometry makes 
working with the L^-norm easier in the Brownian motion case, whereas Poisson processes lend themselves 
naturally to analysis in the L^-norm. 

We remark that it is clear that the choice of scaling laid out in this section and assumed throughout the 
paper will not explicitly cover all cases of interest. For example, one may choose to use approximation 
methods when (i) the abundances of only a strict subset of the constituent species are in an 0{V) scaling 
regime, or [ii) it is the rate constants themselves that are 0{V) while the abundances are 0(1), or (iii) 
there is a mixture of the previous two cases with potentially more than two natural scales in the system. Our 
analysis will not be directly applicable to such cases. However, the purpose of this analysis is not to handle 
every conceivable case. Instead, our purpose is to try and give a more accurate picture of how different 
tau-leaping methods approximate the exact solution, both strongly and weakly, in at least one plausible 
setting and we believe that the analysis detailed in this paper achieves this aim. Further, we believe that error 
analyses conducted under different modeling assumptions can be carried out in similar fashion. 

2.2. Redefining the kinetics. Before proceeding to the analysis, we allow ourselves one change to the 
model detailed in the previous section. As we will be considering approximation methods in which changes 
to the state of the system are determined by Poisson random variables (which can produce arbitrarily large 
values), there will always be a positive probability that solutions will leave a region in which the scaling 
detailed above is valid. Multiple options are available to handle such a situation. One option would be to 
define a stopping time for when the process leaves a predetermined region in which the scaling regime is 
valid and then only perform the analysis up to that stopping time. Another option, and the one we choose, 
is to simply modify the kinetics by multiplying by a cutoff function that makes the intensity functions 
zero outside such a region. This has the added benefit of guaranteeing the existence of all moments of 
the processes involved. Note that without this truncation or some other additional assumption guaranteeing 
the existence of the necessary moments, some of the moment estimates that follow may fail; however, the 



(2.5) 



where < /3 < 1. We note that this scaling satisfies the necessary requirements detailed above as 
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convergence in probability and convergence in distribution results in Theorems 3.10 and 3.17 would still be 
valid. 

Let 7 > be with compact support 17^ C M?J.q, with 7(x) = 1 for all x G Br{x{t)) for some r > 0, 
where x{t) satisfies (2.1). Now we redefine our intensity functions by setting 



(2.6) xl[x)=^[xlV)clJ{-^^^ forxGM^ 

where cjf still satisfies the scaling detailed in the previous section. It is easy to check that the redefined 
kinetics still satisfies (Vx) = VA^{x), where now A^{x) has also been redefined by multiplication by 
7(x). Further, the redefined is identical to the previous function on the domain of interest to us. That 
is, they only differ if the process leaves the scaling regime of interest. For the remainder of the paper we 
assume our intensity functions are given by (2.6). Finally, we note that for each k we have the existence of 
an Lfc > such that 

(2.7) sup < Lk. 

xi^W-, \a\<oo 

3. Strong error analysis for Euler and midpoint tau-leaping. Throughout this section we assume a 
time discretization = Iq < ti < ■ ■ ■ < t]\i = T with t„ — = = for some < /? < 1. In 
Section 3. 1 we give some necessary technical results. In Section 3.2 we give bounds for supj<y E|X^(t) — 
Z^{t)\ and supi<rE|X^(t) - Z^{t)\ in terms of V, where X^{t), Z^{t), and (t) are the normalized 
processes and satisfy the representations 

(3.1) = X^(0) + ^ E^*^ [v j\l{X^{s))ds^ uk 

(3.2) = X^(0) + :i ^ Yu [v 1^ AliZ"" o r^[s))ds^ 

(3.3) Z^it) = X^(0) + :^ n [v 1^ AlopVoZ^o r]{s)ds^ Uk, 



where 



2 ^--k[z)iyk, 

k 



and r]{s) = tn for s E tn+i). In Sections 3.3 and 3.4 we use different couplings of the processes than 
those above to provide the exact asymptotics of the error processes — and X^ — Z^ . 

3.1. Preliminaries. We present some technical, preliminary concepts that will be used ubiquitously 
throughout the section. For a more thorough reference of the material presented here, see [6], chapter 6. 
We begin by defining the following filtrations that are generated by the Poisson processes 1^, 

Tu = a{Yk{sk) : Sk < Uk} 

Tl, = a{Yk{r),Yds) : /c / i, s < n, r G [0, oo)}, 
where n is a multi-index and u is a scalar. 
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Lemma 3.1. Suppose that X{t) satisfies (1.2) with non-negative intensity fiinctions A^. For t > and 
a choice of k, 

(3.4) Tfc(t)= [\k{X{s))ds 







is an {T^} -stopping time. 

Proof. For u > let a{u) satisfy 

fa{u) 



\k{X{s))ds = u, 

where we take a{u) = oo if Xk{X{s))ds < u. Then a{u) is adapted to and {rfc(t) < u} = {t < 
a{u)] e Ft □ 

Therefore, if the processes X{t) and Z{t) satisfy (1.2) with non-negative intensity functions A^^i and 
Afc^2> respectively, then for t, s > and a choice of k 



(3.5) 

E 



Yk[ 1^ Xk,i{X{r))drj-Yki^J^ \k,2{Z{r))dr 



= E 


L 







ft 



Xk,i{X{r))dr - / \k,2{Z{r))dr 







because (i) both the maximum and minimum of two stopping times are stopping times, and {ii) Yk is 
monotone. 

Similarly to above, one can show that T{t) = (Ti(f), r2(t), . . . ), where rfc(t) is as in (3.4), is a multi- 
parameter {J^ti} -stopping time. We now define the filtration 

yt = >T(t) 

and note that by the conditions of Section 2.2 the centered process 

t \ / ft \ ft 

def 



(3.6) Yk i^J^ Xk{X{s))dsj =YkiyJ^ XkiX{s))dsj - A,.(X(s))(is, 

is a square integrable martingale, with respect to Qt, with quadratic variation (^J Xk{X{s))ds^ . This 
fact will be used repeatedly throughout the paper. 

3.2. Bounds on the strong error. The following theorems give bounds on the errors supj<2^ E|X^(t) — 
andsupi<7.E|X^(f) 

Theorem 3.2. Let X^ {t) and {t) satisfy (3.1) and (3.2), respectively, for t < T. Then there exists 
a constant C = C{T) > such that 

supE|X^(t) - Z^{t)\ < CV-^ = Ch^. 

t<T 
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Proof. Forte [0, T] define = E|X^ (t) - Using (3.5) and (2.7) 

E{t) < (V Wk\Lk)E f \X^{s) - o r,{s)\ ds 
k 

<{y2Wk\Lk) [ E{s)ds + {y2\uk\Lk)E [ \Z^{s)-Z^ orj{s)\ds. 
k k 

Tlie second term on tlie rigiit above can be bounded similarly 

E f \Z^{s) - Z^ o rj{s)\ ds < (V luklLkW-f", 

k 

and the result holds via Gronwall's inequality. □ 

Theorem 3.3. Let {t) and {t) satisfy (3.1) and (3.3), respectively, for t < T. Then there exists 
a constant C = C{T) > such that 

supE|X^(t) - Z^{t)\ < CV-<'^\ where k{(3) = min|i^,2/3j . 
t<T I 2 J 



Before proving Theorem 3.3 we present some preliminary material. Let = and 

define 

U^^^{s) = Z^{s) -p^ oZ^ o r,{s) = Z^{s) -Z^o rj{s) - ]^V-^F^{Z^ o r]{s)) 

and 

(s - ry(s) - o ry(s)). 

Then 

(3.7) [/^■^(s) - U^^^{s) = Z^{s) -Zy o r,{s) + {s - rj{s)){F^{p^ o o 7]{s)) - {Z^ o r,{s))), 
where Z^{t) = Z^{t) - F^{p^ o Z^ o r]{s))ds is a martingale. 
Lemma 3.4. For all < 13 < 1, there exists a C > such that 

sup E|C/^'i(s) - U^'\s)\ < C7y-'=(^). 

s<oo 



Proof. Clearly, the third term on the right of (3.7) is 0{V '^^) uniformly in s. Thus, 

E\U^^\s) - U^'\s)\ < E\Z^{s) -zy o T]{s)\ + 



< 



-Vlz^fcl^E/ Al{p^ oZ^ or]{T))dT\ +ciV 



1/2 

2/3 



/r?(s) 

<C2y-(l+/^)/2 + ciy-2/^, 

for constants ci and C2 which do not depend upon s. □ 
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Lemma 3.5. For alio < (3 <1 and < t, and for a E {2, 3, 4, ... } 

lim supE[\U^'\s) - U^'HsT] = 0. 

Proof. The third term on the right of (3.7) is 0{V~'^l^), so 



tVA 



-4/3 ^ 



< 



V 



Izvfcl^E / Al{p^ oZ^ o r]{r))dr + CV' 



0{V 



-{(l+/3)A4/3)^ 



showing the a = 2 case. 

It is simple to show that sup^<j E[|f7^'^(s) — f/^'^(s)|"] is uniformly bounded in V for any a G Z>o. 
The a = 2 case then gives the necessary bounds for the arbitrary a case. □ 

Note that by Lemmas 3.4 and 3.5 

Al{Z^{s)) - Alip'' o o r^{s)) = VAr(/ o o • ^/^'^(s) + ©(^-^Z^) 
(3.8) = VAlip"^ o o r]{s)) ■ U^^^{s) + 

We finally note that for any bounded function g and any n > 



5(r?(s))C/^'^(s)(is = 



and so for any t > 



(3.9) 



g{'n{s))U''^\s)ds 



{2t - 2r]{t) - V 



V 



-2/3 



5(r/(t))F^(Z^o^(t)) = 0(y-2/3) 



Proof. (Of Theorem 3.3) For t < T define E{t) = E\X^ (t) - {t)\. Letting a denote constants 

E{t)<Y,Wk\E f Al{X''{s))ds- f Alop^iZ"" or^{s))ds 
k Jo 

<ci f E{s)ds + y^\uk\E f Al{Z^{s))-AXop^{Z'' or^{s))ds 
Jo f, Jo 

<ci [ E{s)ds + C2V-^'^'^\ 
Jo 

where the final inequality used both (3.8) and (3.9). The result now follows from Gronwall's inequality. □ 

3.3. Exact asymptotics for Euler tau-leaping. Throughout this section and the next all convergences are 
understood to hold on bounded intervals. More explicitly, we write — )• X if \\mv-^oo P{s^Vt<T I (*) ~ 
> e} = for all e > and T > 0. Because of the simplifying assumptions made on the kinetics in 
Section 2.2 it is not difficult to show that X^ — X also implies limv^oo Esup^<7i \ X'^ {t) — X{t)\ =0. 
In light of this, when we write X^ = + 0{V~^) for some p > in this section and the next we mean 
that for any T > there exists a C{T) such that 



lim FPEsup|X''(t) - Z''(t)| < C7(r). 
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Finally, recall that F^{x) = Efc^fe (^)^fc 

and note that the function F{x) and the deterministic process 
x{s) used in the characterization of the error processes are defined via (2.1). 

Theorem 3.2 suggests that — scales like V~^. In this section we make this precise by characteriz- 
ing the limiting behavior of V^{X'^ — Z^), as F — )• oo. To get the exact asymptotics for the Euler tau-leap 
method we will use the following coupling of the processes involved 



(3.10) 



V/r^V 



+ Yk,2 ( V I AliX^'is)) - A', {X' (.)) A Al {Z' o r,{s))ds 



V tvV I 



\V (ryV 



(3.11) 



YkAV I Al{X''{s))AAl{Z' o^{s))ds 



V ( r,V 



+ n,3 [VI AliZ'^o r^is)) - Al iX' (s)) A A^ {Z' o r?(s))ds 



V tvV I 



\VlryV 



It is important to note that the distributions of X^ and Z^ defined via (3.10) and (3.11) are the same as 
those for the processes defined via (3.1) and (3.2). 

The following Lemma is easy to prove using Doob's inequality. 



Lemma 3.6. For and given by (3.10) and (3.11), 



7V 



0. 



Combining Lemma 3.6 and (2.2) shows that Z^ 



0, where x is the solution to the associated ODE. 



Similarly Z^ o t] — x ^ 0. These facts will be used throughout this section. 
Centering the Poisson processes, we have 



X^{t) - Z^{t) = M^{t) + f F^{X^{s)) - F^{Z^ o r]{s))ds 

Jo 

= M^{t)+ [ F^ {X^ is)) - F^ {Z^ {s))ds + [ F^{Z^{s))-F^{Z^ or]{s))ds, 
Jo Jo 



(3.12) 



where is a martingale. 

To obtain the desired results, we must understand the behavior of the first and third terms on the right of 
(3.12). We begin by considering the third term. We begin by defining C/^'^ and U^''^ by 



U^'\s) = - o r]is), U^''{s) = {s- r]{s))F^ {Z" o r/(s)). 



7V 



V,2/ \ drf 



(ryV 



Then, 



Uy^\s) - U'^'^s) = Z^'is) - Zy o ry(s), 
where = Z^ (i) - o r?(s))ds is a martingale. Thus 

F^{Z^{s)) - F^{Z^ o ry(s)) = DF^{Z^ o i^{s))U^^'^ {s) + 0(^-^/3) 

(3.13) = DF^{Z^ o r/(s))?7^'2(s) + DF^{Z^ o rj{s)){U^'\s) - U^'^{s)) + 0{V- 

Lemma 3.7. For a// < /3 < 1, < t, and a G {2, 3, 4, ... } 



lim y"^supE[|C/^'^(s) - C/''''(s)n = 0. 
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Proof. The proof is similar to that of Lemma 3.5. □ 
We may now characterize the limiting behavior of the third term of (3.12). 
Lemma 3.8. For < /3 < 1 and any t > 0, 

yP f F^{Z^{s))- F^{Z^ oi^{s))ds^- f DF{x{s))F{x{s))ds. 
Jo 2 Jo 

Proof. By (3.13) and Lemma 3.7 

/* F^{Z^{s)) - F^{Z^ o r^{s))ds = f DF^{Z^ o r]{s))F^{Z^ o r]{s)){s - 7]{s))ds + eY{t), 
Jo Jo 

where e]^ — )• as 1/ — )• oo. By Lemma 3.6 convergence results similar to (2.2) hold for the process Z^ o r], 
and because f^l^g^~^^ {r — rj{s))dr = ^V~'^^, the lemma holds as stated. □ 

Turning now to AI^ , we observe that the quadratic covariation is 



y2 

k 



where 



Nlsit) = Yk (v f AliZ"" o r^{s)) - A AI{Z'' o r^{s)) 







which as y — 5- oo is asymptotic to 

(3.14) i ^ r \Al{X''{s)) - AliZ"" o r^{s))\ds 

u •'0 



k 

We have the following Lemma. 

Lemma 3.9. For < /3 < 1, V^M"^ asV ^ oo. 

Proof. Multiplying (3.12) by F", we see that V'^{X^ - Z^) provided a < /? (so that the third 
term on the right goes to zero) and provided V^M^ — )• 0. By the martingale central limit theorem, the 
latter convergence holds provided y^°[M^] — (see Lemma A.2 in Appendix A). Let ao = sup{a : a < 
P, F^"[M^] — > 0}. Since ao < 1, we have that 2ao — 1 < cuo < /3, which implies by the definition of ao 
that y2ao-i(x^ - Z^) 0. Therefore 

''ki^k 



V^""[M%^y] [ V^''^-^\Al {X^ (s)) - Al {Z^ o 7]{s))\ds 
k -^0 

^ E r ^'""^'iV^r (^'' o vis)) ■ (Z'^is) - Z^ o r^{s))\ds VkV, 
k 

« V /" y2ao-i|vA^(Z^ o 7?(s)) • F^{Z^ o r]{s))\{s - r]{s))ds 
Jo 
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where in the second approximation we used that y'^<^o-i(^x^ — Z^) — )• 0, in the third approximation we 
substituted f/^'^(s) for [/^'^(s), and by / gf we mean / — s-Oasy— s-oo. The last expression goes to 
zero whenever 2aQ — 1 < /?, hence the convergence holds. □ 

We now have the following theorem characterizing the behavior of V^{X^ — Z^\ 

Theorem 3.10. ForX"^ and Z^ given hyO-^^) and O.W) and for ^ < ^ < 1, V^{X^ -Z^) ^ £, 
where E is the solution to 

(3.15) £{t)= I DF{x{s))£{s)ds + \ I DF{x{s))F{x{s))ds, £{0) = 0. 

Jo 2 Jo 



Proof. Multiply (3.12) by ^ and observe that 

f F^ {X^ [s)) - F^ {Z^ {s))ds ^ f DF^{Z^{s))V^{X^{s)-Z^{s))ds. 
Jo Jo 

The theorem now follows directly from Lemmas 3.8 and 3.9. 



□ 



3.4. Exact asymptotics for midpoint tau-leaping. Throughout this section the Hessian matrix associated 
with a real valued function g will be denoted by Hg. Also, for any vector U, we will denote by U^HF^ {^W 
the vector whose ith component is HFYU, and similarly for F. 

The goal of this section is to characterize the limiting behavior of V^^^^X^ {t) — Z^{t)) where 



k{I3) = min{2/3, (1 + /3)/2} 



2/3 /? < 1/3 

(l + /3)/2 /3>l/3 • 



To get the exact asymptotics for the midpoint method we will use the following representation of the pro- 
cesses involved. 



(3.16) X^(t) = X^(0) + i5] ^"^'{^l 



n,i ( V I AliX^'is)) A Al (p^ o o rj{s)) ds 



+ n,2 [V AliX^'is)) - AliX^'is)) A o Z"" o ^{s))ds 

(3.17) Z^(t) = X^(0) + i^ Y,,, (v I^UX'' (s)) A Alip'' o Z"" o ^{s)) ds^ 

+ n,3 (v Alip^ o Z^ o r^{s)) - AliX^'is)) A Al{p^ o Z^ o r^{s))ds 
The following is similar to Lemma 3.6. 

Lemma 3.11. For X^ and Z^ given by (3.16) and (3. 17), X^ - Z"^ ^ 0. 

Combining Lemma 3.11 and (2.2) shows that Z^ — x — 0, where x is the solution to the associated 
ODE. Similarly Z^ o ry — x — )• 0. These facts will be used throughout this section. 
Centering the Poisson processes, we have 

X^{t) - Z^{t) = M^{t) + f F^{X^{s)) - F^ip^ oZ^ o rj{s))ds 

Jo 

(3.18) =M^{t)+ [ F^ {X^ (s)) - F^ {Z^ {s))ds + [ F^ {Z^ (s)) - F^ {p^ o Z^ o r]{s))ds, 

Jo Jo 
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where is a martingale. 

As before, we must understand the behavior of the first and third terms on the right of (3.18). We begin 
by considering the third term. Proceeding as in the previous sections, we define U^'^ and U^'^ as 

= Z^{s) -p^ oZ^ o rj{s) = Z^{s) -Z^ o r]{s) - U^-^F^{Z^ o r]{s)) 

and 

(jV^^i^s) ^' {s - r^is) - o r?(.)). 

Then 

(3.19) U^'^is) - U^'^is) = Z^{s) -Z^o r]{s) + {s - r]{s)){F^{p^ o Z^ o n{s)) - {Z^ o rj{s))), 
where Z^ (t) = Z^ (t) - F^ {p^ o Z^ o r]{s))ds is a martingale. Then 

F^{Z^{s)) - F^{p^ o Z^ o 7]{s)) = DF^{p^ o Z^ o r]{s))U^'^{s) 

+ ]^U'''\sYhf''{p^ o o ii{s))U''^\s) + 0{V~^^) 
= DF^ {p^ o Z^ o n{s))U^'^{s) 

(3.20) + ^U^'^{sfHF^{p^ oZ^ o r]{s))U^'^{s) 

+ DF^{p^ oZ^ o rj{s)){U^'^{s) - U^'^is)) + 0{V'^l^). 

Lemma 3.12. For all < /3 < 1, < t and a G {2, 3, 4, ... } 

lim y"'^supE[|C/^'3(s) - U^'^{sT] = 0. 

V-S>oo 

Proof. The proof is similar to Lemma 3.5. □ 
Let 

KiiP) = min{2/3,/3 + l/2} = 
Note that > k{P) for all (3 > 0. 

Lemma 3.13. For < /? < ^ and each t > 0, 

(3.21) f' DF^{p^ oZ^ or]{s)){U^'^{s)-U^'^{s))ds 
Jo 

for (3 = I 

(3.22) V f' DF^{p^ o Z^ or]{s)){U^'^{s) - U^'^{s))ds Mi{t) + - f DF{x{s)fF{x{s))ds, 

Jo 4 Jo 

and for ^ < /? < 1, 

(3.23) f' DpV^^pV ^ ^ j^(s))([/V,3(g) _ [7V,3(^))^^ ^ 

Jo 

where Mi is a mean zero Gaussian process with independent increments and quadratic covariation 

(3.24) [Mi]t = ^ l''Y,Ak{x{s))DF{x{s))ukU^DF{xis)fds. 

•^0 k 
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t 

DF{x{s)fF{x{s))ds, 
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Proof. By Lemma A. 1 in Appendix A, 

MY{t) = f DF^ip^ oZ^o 7]{s)){Z^{s) -Z^ o 7]is))ds 
Jo 

+ DF^ip^ o Z^ o r]{t)){Z^{t) - Zy o 7?(t))(r?(t) + F"^ - £) 
is a martingale and its quadratic covariation matrix is 

\\r]{s) + V-^ - sfDF^{p^ oZ^ o r]{s))d[Z^]sDF^{p^ o Z^ o i^{s)f . 
Jo 

Noting that J^l^^)^^'" {r]{r) + F"^ - rfdr = ^V-^l^, it follows that 

V^^+^[MY]t^\ fTA,ix{s))DF{x{s)Hi^'[DF{x{s)fds, 

so by the martingale central limit theorem V^'^^^'^M^ converges in distribution to a mean zero Gaussian 
process with independent increments and quadratic variation (3.24). 

Since F V2(^v _ o r]) ^ 0, the integral on the left of (3.21), (3.22) and (3.23) can be replaced by 

(3.25) MY{t) + r DF^ip^ oZ^o rj{s)){s - i]{s)){F^{p^ o Z^ o ry(s)) - F^ {Z^ o r,{s)))ds 
Jo 

without changing the limits. The second term in (3.25) multiplied by V^^ converges to ^ DF{x{s))'^F{x{s))ds 
on bounded time intervals and the three limits follow. □ 

Lemma 3.14. For < /? < 1, 

y2/3l r u^,3{s fHF^{p^ oZ^ or]{s))U^'\s)ds^ ^ [' F{x{s)f HF{x{s))F{x{s))ds. 

Proof. By Lemma 3.12 we canreplacef/^-^byC/^'^. Observing that f^l^^^^^ \s-r]{s)-lV~'^)'^ds = 

12 ^ 

y'^^\ [\s - V(.s) - \v-^fF^{Z^ o 7]{s)fHF^{p^ oZ^o 7^{s))F^{Z^ o r?(s))ds 
2 Jo 2 

converges as claimed. □ 
We may now characterize the behavior of the third term of (3.18). 
Lemma 3.15. Let 

R^{t)= / {s-r,{s)- -V-^)DF^{p^ oZ^ on{s))F^{Zor]{s))ds. 
Jo 2 

Then for < /? < 

y2/3 (^J^ {F^{Z^{s)) - F^{p^ o Z^ o r]{s))) ds - R^{t)j 

^ ^* DF{x{s)fF{x{s))ds + ^J^ F{x{s)fHF{x{s))F{x{s))ds, 
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16 ANDERSON, GANGULY, AND KURTZ 

for^ = \, 

Mi{t) + \ I DF{x{s)fF{x{s))ds + 7^ / F{x{s)YHF{x{s))F{x{s))ds 
and for ^ < /3 < 1, 

yP+i/2 f\Fy[z^{s)) - F^{p^ oZ^o rj{s)))ds Mi{t). 
Jo 

Remark. Note that V^'^R^h uniformly bounded, or] = 0, and 

R'^it) = l[it- V{t)f -{t- mV-^] DF^'ip'' o o r/(t))F^(Z o r?(t)). 

Proof. The lemma follows from (3.20), the previous lemmas, and by noting that Jq DF^ {p^ o Z^ o 

ri{s))U^'^s)ds = R^{t). □ 

We now turn to and observe that 

k 

where 

<2(i) = Yk (v 1^ AliX^'is)) - AliX^'is)) A AUp'' o o r^{s))ds 
NU^) = Yk (v 1^ Alip'' o o ^{s)) - AliX'^is)) A o o ^{s))ds 

which as y — )• oo is asymptotic to 

Consequently, we have the following. 



A ■ 



Lemma 3.16. for < /3 < 1, f(i+^)/2m^ ^ M w/jere M is a mean-zero Gaussian process with 
independent increments and quadratic covariation 



[^]* = E 4 |VAfc(x(s)) • F{x{s))\ds Uk^l. 



k 

Proof. Multiplying (3.18) by ^ , we see that V^iX^ -Z^) provided a < (so that the third 

term on the right goes to zero) and provided V^M^ — 0. By the martingale central limit theorem, the latter 
convergence holds provided V'^°'[M^] 0. Let ao = sup{a : a < (/3 + l)/2, V'^"[M^] 0}. We make 
two observations. First, because oq < 1, we have that 2ao — 1 < oq. Second, because ao < (/? + 1) /2, we 
have that 2ao — 1 < /3, and, in particular, 2qo — 1 < ki{P) for all (3 G (0, 1). Combining these observations 
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with the definition of shows that y2(2Qo-i)[j^y]^ _^ q ^nd hence y2ao-i(-x^ - Z^) 0. We now 
have 



T 



V'"^'[M%^y] [ V^^^^~'\Ak{X''is))-Akip''oZ''orj{s))\ds 
k 

^y^f y2ao-i|, _ ^(,) _ ly-/3||vAfc(p^ o o ^(s)) ■ o r,{s))\ds u^ul 

k ^ 

where in the second hne we used that l/^ao-i^^v _ _^ ^j^^ jj^^j^ substituted [/^'^(s) for f7^'^(s). 
Since the last expression would go to zero if 2ao — 1 were less than /?, we see that 2aQ — 1 = (3, that is, 
«o = (/? + l)/2- Furthermore, observing that /J^^j)^"''^ |s — '^(•s) — = \V~'^^ , we see that 

k •'^ 

and the lemma follows by the martingale central limit theorem. □ 
Collecting the results, we have the following theorem. 
Theorem 3.17. Let 

n{t) = \ [ DF{x{s)fF{x{s))ds + ^ [ F{x{s)fHF{x{s))F{x{s))ds 
ForO< P <l, -Z^ - R^) £i, where 8i is the solution of 

(3.26) £i{t)= I DF{x{s))£i{s)ds + n{t), fi(0) = 0. 

Jo 

For 13 = \, -Z"^ - R^) £2, where £2 is the solution of 

(3.27) £2{t) = M{t) + [ DF{x{s))£2{s)ds + n{t), £2(0) = 0. 

Jo 

Forl< 13 <1, y{i+/3)/2(x^ - Z^) £3, where £3 is the solution of 

(3.28) £3it) = M{t) + / DF{x{s))£3{s)ds, £3(0) = 0. 

Jo 



Proof. For P < I, R^ is 0{V~'^^). Subtract R^ from both sides of (3.18) and observe that 

r F^{X^{s)) - F^{Z^{s))ds « f DF^{Z^{s)){X^{s) - Z^ (s) - R^ {s))ds 
Jo Jo 

+ / DF^{Z^{s))R^{s)ds. 
Jo 

y2P f' jjpV^^V^^y^j^V^^^^^^ _}_^ [' DF{x{s)fF{x{s))ds, 
Jo 12 Jo 

imsart-aap ver. 2009/02/27 file: AndGangKurtzRevision.tex date: July 9, 2010 



Since 



18 



ANDERSON, GANGULY, AND KURTZ 



the first two parts follow from Lemmas 3.15 and 3.16. 

For /3 > ^, (1 + /3)/2 < 2/3 A so y(i+/3)/272V _^ q 

y(i+p)/2 f' F^'iZ^^s)) - F^ip"" o o ^{s))ds ^ 0, 
Jo 

and the third part follows by Lemma 3.16. □ 

4. Weak error analysis. As in previous sections, we assume the existence of a time discretization 
^ = ti) < ti < ■ ■ ■ < ti^ = T with tn - tn-i = for some < /? < 1. We also recall that ri{s) = tn for 
in < s < tn+i for each n < N — I. 

Let be a Markov process with generator 

(4.1) iA''f){x) = Y,VAl{x)if{x + Uk/V) - fix)). 

k 

Defining the operator 

(4.2) {B^J){x) = J2 V^ki^Kfix + ^k/V) - fix)), 

k 

we suppose that and are processes that satisfy 

(4.3) = E/(Z o r?(t)) + E T 

Jr,{t) 

and 

(4.4) = EfiZ o ry(t)) + E T (Vo2or,(t)/)(^(5))d^, 

J»,{t) 

for all t > 0, respectively. 

We begin with the weak error analysis of Euler tau-leaping, which is immediate in light of Theorem 3.10. 

Theorem 4.1. Let it) be a Markov process with generator (4.1) and let it) be a process that 
satisfies (4.3) for the operator (4.2). Then, for any continuously differentiable function f and any t <T 

lim (E/(X^(t)) - E/(Z^(t))) = Sit) ■ Vfixit)), 



where Sit) satisfies (3.15). 

Proof. Without loss of generaUty, we may assume that X^it) and Z^it) satisfy (3.10) and (3.11), 
respectively. The proof now follows immediately from a combination of Taylor's theorem and Theorem 
3.10. □ 

Remark. Because the convergence in Theorem 4.1 is to a constant independent of the step-size of the 
method, we see that Richardson extrapolation techniques can be carried out. However, we have not given 
bounds on the next order correction, and so can not say how much more accurate such techniques would be. 

We now consider the weak error analysis of the midpoint method. 
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Theorem 4.2. Let [t) be a Markov process with generator (4.1) and let {t) be a process that 
satisfies (4.4) for the operator (4.2). Then, for any two times continuously differentiable function f with 
compact support, there exists a constant C = C{f,T) > such that 

V'^^\Ef{X^{T)) - Ef{Z^{T))\ < C. 



Before proving Theorem 4.2, some preliminary material is needed. Let = {y '■ y = x/V, x £ Z'^}, 
and for x S LX and a given function /, let 

(4.5) v{t,x)=E^f{X^{t)), 

where represents the expectation conditioned upon X^{0) = x. Standard results give that v{t, x) satisfies 
the following initial value problem (see, for example, [7] Proposition 1.5) 

(4.6) dtv{t,x) = A^v{t,x) = ^VA'^{x){v{t,x + i^k/V) -v{t,x)), v{0,x) = f{x), xel7. 

k 

The above equation can be viewed as a linear system by letting x enumerate over and treating v{t, x) = 
Vx{t) as functions in time only. It can even be viewed as finite dimensional because of the conditions on the 
intensity functions A^. That is, recall that Af^{x) = for all x outside the bounded set (see Section 
2.2); thus, for any such x ^ v{t, x) = Vx{t) = f{x), for all t > 0. 

For concreteness, we now let M denote the number of reactions for the system under consideration. For 
k,ee [1,...,M] andrc G L^let 

(4.7) Dk{t, x) = V{v{t, X + Uk/V) - v{t, x)) 

(4.8) Dkeit, x) = V{Dk{t, x + u,/V) - Dk{t, x)). 

represent approximations to the first and second spatial derivatives of v{t, x), respectively. For notational 
ease, we have chosen not to explicitly note the V dependence of the functions v{t, x), Dk{t, x), or Df^£{t, x). 

The following lemma, which should be viewed as giving regularity conditions for v{t, x) in the x variable, 
is instrumental in the proof of Theorem 4.2. The proof is delayed until the end of the section. 

Lemma 4.3. Let v{t, x), Dk{t, x), and Dke{t, x) be given by (4.5), (4.7), and (4.8), respectively, and 
let T > 0. There exists Ki > and K2 > that do not depend upon V such that 

(4.9) sup sup sup \Dk{t,x)\ < Ki 

t<T k<M xGhV 

(4.10) sup sup sup \Dki{t,x)\ < K2. 

t<T k/<M x&LV 



We will also need the following lemma, which gives regularity conditions for Dk{t, x) in the t variable, 
and whose proof is also delayed. 

Lemma 4.4. Let Dk{t, x) be given by (4.7). There exists a K > that does not depend upon V such 
that 

sup sup sup \Dk{t + h,x) — Dk{t, x)\ < Kh, 

t<Tk<Mx&y 

for all h > 0. 
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Proof, (of Theorem 4.2). Define the function u{t, x) : [0, T] x R by 

(4.11) u{t,x)=EJ{X^(T-t)), 

and for any w{t, define the operator £ by 

Cw{t,x) = dtw{t,x) +A^w{t,x) = dtw{t,x) + '^VAl{x){w{t,x + Uk/V) - w{t,x)). 

k 

Note that u{t, x) = v{T — t, x), where v{t, x) is given by (4.5), and so by (4.6) Cu{t, x) = for t G [0, T] 
and X € L^. We also define the operator 

Czw{t, x) = dtw{t, x) + B^w{t, x) = dtw{t, x) + ^ V A)^{z){w{t, x + fk/V) — w{t, x)), 

k 

so that by virtue of equation (4.4), for t < T and any differentiable (in t) function w{t, x) 



(4.12) Ew{t, Z^{t)) = Ew{r,{t),Z^ o 7]{t)) + / ECpV,zVorj(t)w{s, Z^{s))ds. 

Jr){t) 

Recalling (4.11), we see that 



Eu{T, Z^{T)) = Ef{Z^{T)) 

Eu{T,X^{T)) = Eu{0,X^{0)) = Ef{X^{T)). 

Therefore by (4.12), and using that X^{0) = Z^{0), 

Ef{Z^{T)) - Ef{X^{T)) = Eu{T, Z^ (T)) - Eu{0, Z^ {0)) 

N-l 



Eu{tn+i,Z^{tn+i)) - Eu{tn, Z^itn)) 

n=0 

CpV^zv(^t^~,u{s,Z^{s))ds. 

n=0 



Because Cu{t, x) = for t < T and x £ L 



V 



rtn+i rtn+1 

^ ^pVoZV{t„)uis,Z^{s))ds = E CpVozv(t„)uis,Z^{s)) - Cu{s,Z'^{s))ds 

= Ye V (p^ o - AliZ^'is))] {uis, Z^'is) + ,^k/V) - uis, Z^is))) ds 

k 

(4.13) = J^E [A]:ip''oZ''itn))-AliZ^is))]DkiT-s,Z''{s))ds. 

1.. J tr) 



Thus, it is sufficient to prove that each of the integrals in (4. 13) are 0{V '^1^). By Lemma 4.4 each integral 
term in (4.13) can be replaced by 

(4.14) ll{trr)=E / [Al{p^ oZ^'itn))- Al{Z''{s))\Dk{T -tn,Z^{s))ds. 

J tn 

The remainder of the proof consists of proving that l]^ (tn) = 0{V^^'^). 
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Letting g}^{x) = [A^ip^ o (tn)) - A^ix)] Dk(T - t„,x) and applying (4.4) to the integrand in 
(4.14) yields 

rtn+i 



rtn+1 

I^iU) =E [yl^/ o Z^{tn)) - AliZ^'itn))] Du{T - U, {tn))ds 

J tn 

+ Y.^ / ^^(p^ o Z^(t„))(5r(2^(r) + y,/V) - g"" {Z"" {r)))drds. 



We have ^^(p^ o = + V A^ {Z'' (1,,)) • '^V'^ {Z'' {tn))u, + 0{V-^P). 

Thus, 



(4.15) 



1 /■tn+1 

(4.16) +E^/ / VAj{p''oZ''{tn)){gX{Z''ir) + u,/V)-g''{Z^{r)))drds 



3/3^ 



After some manipulation, the expected value term of (4.16) becomes 

E / VAj{p''oZ''{tn))[Al{Z''{r))-Al{Z''{r) + u,/V)]Dk{T-t^,Z''{r) + u,/V)drds 

J tji J t fi 

+ E / Aj{p''oZ''itn))[Al{p''oZ''itn))-AliZ''{r))]Dk,{T-tn,Z''{r))drds. 

'J tn 'J tn 

By Lemma 4.3 the last term above is 0{V~'^l^). Taylor's theorem and the fact that A^ {p^^ o Z^ {tn)) = 
Ay {Z^ {tn) + 0{V~^) then shows us that the expected value term of (4.16) is equal to 



-E / A]{Z^{tn))VAl{Z^{r)) • u,Du{T - tn, Z"" {r) + yj/V)dr ds + 0{V-''^) 

J tn ^ tn 

rtn+1 rs 

(4.17) =-E/ / Aj{Z^{tn))VAl{Z^{r))-i,,Dk{T-tn,Z''{r))drds + 0{V-'"'), 

where the second equality stems from an application of Lemma 4.3. 

By Lemma 4.3, the function = Aj {Z^ {tn))VA^ {x)-UjDk{T-tn, x) satisfies sup^ \(j){x+ue/V)- 
4>{x)\ = 0{V~'^). Therefore, applying (4.4) to (4.17) shows that (4.17) is equal to 



-E 



/ / Aj{Z^{tn))VAl{Z^{tn))-i^,Dk{T-tn,Z^{tn))drds + 0{V-^^) 

J tn ^ tn 



Noting that the sum over j of the above is the negative of (4.15) plus an 0{V ^^) correction concludes the 
proof. □ 

Theorem 4.2 can be strengthened in the case of /3 < 1/3. 

Theorem 4.5. Let {t) be a process with generator (4.1) and let Z^ {t) be a process that satisfies 
{A A) for the operator (4.2). Suppose also that (3 < 1/3. Then, for any continuously differentiable function 
f, 

lim (E/(x^(r)) - E/(^^(r))) = £^{t) ■ v/(x(r)), 

where Si{t) satisfies (3.26). 
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Proof. Noting that (T) = 0, this is an immediate consequence of Theorem 3.17. □ 

Remark. In Theorem 4.1 we provided an explicit asymptotic value for the scaled error of Euler tau- 
leaping in terms of a solution to a differential equation for all scales, < /3 < 1, of the leap step. However, 
Theorem 4.5 gives a similar result for the midpoint method only in the case < /3 < 1/3. For the case 
1/3 < /3 < 1, Theorem 4.2 only shows that the error is asymptotically bounded by a constant. The reason 
for the discrepancy in results is because in Section 3 we were able to show that the dominant component of 
the pathwise error for Euler tau-leaping for all /3 G (0, 1) and for midpoint tau-leaping for /3 € (0, 1/3) was 
a term that converged to a deterministic process. However, in the case /3 > 1/3 for midpoint tau-leaping, the 
dominant term of the error is a non-zero Gaussian process. We note that this random error process should 
not be viewed as "extra fluctuations," as they are present in the other cases. In these other cases, they are just 
dominated by the error that arises from the deterministic "drift" or "bias" of the error process. We leave the 
exact characterization of the weak error of the midpoint method in the case /3 > 1/3 as an open problem. 

We now present the delayed proofs of Lemma 4.3 and Lemma 4.4. 

Proof. (Of Lemma 4.3) Let Ci > be such that 

sup sup |i:>fc(0,x)| = sup sup \V{f{x + i^k/y) - fix))\ < Ci. 

Using (4.6), a tedious reordering of terms shows that Dk{t, x) satisfies 

dtDk{t, x) = Y^ Aj{x)V[Dk{t, X + Uj/V) - Dk{t, x)] 

(4.18) 

+ Y.^A){x + Vk/V) - A]{x))VDj{t,x + Uk/V). 

3 

Similarly to viewing f (t, x) = v^if) as a finite dimensional hnear system, (4.18) can be viewed as a linear 
system for the variables -Dfc(t, x) = -Dj^ ^.^(t), for /c S [1, • • • , M] and x E L^. Because Aj(x) = for all 
X ^ Vl^, we see that dtDj^{t, x) = for all x such that x ^ Q,j and x + i^j/V ^ for all j G [1, • • • , M]. 
Therefore, the system (4.18) can be viewed as finite dimensional also. 

Let Li = [1, . . . , M] X L^. We enumerate the system (4. 18) over 6 G Ti. That is, for b = {k, x} G Ti we 
let -Db(t) = Dk{t, x) = Dh^ [t, 62)- After some ordering of the set Fi, we let MF^ denote the set of (infinite) 
vectors, v, whose b^^ component is Vb G M, and then denote D{t) G M^^ as the vector whose 6*^ component 
is Db{t). Next, for each b = {k, x} G Fi, we let 

3 3 

and let Rf, G satisfy 

Rb ■ V = ^Aj{b2)v{b,,b2+uj/v} 
j 

n-V = Y,{A]{b2 + Vbjy) - A]{b2))Vv^^^b2^,^Jy^^, 
3 

for all V G M'"^ . It is readily seen that for any b both Rb and Vb have at most M non-zero components. 
Also, by the regularity conditions on the functions ^]^'s, the absolute value of the nonzero terms of are 
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uniformly bounded above by some K, which is independent of V. Finally, note that Ri,-1 = S^. Combining 
the previous few sentences shows that for any vector v G MF^ we have the two inequalities 



(4.19) \Rb-v\ 



j 

(4.20) \n ■ v\ < KM\\v 



< Sb\\v\ 



I OO ) 







where ||u||oo = sup^gpi I'^bl- We now write (4.18) as 

D'^,{t) = -VSbDbit) + VRb ■ D{t) + n ■ D{t), 

and so for each 6 G Fi 

(4.21) ^Dbitf = -2VSMtf + 2VDb{t)Rb • D{t) + 2Db{t)rb • D{t). 

at 

Only a finite number of the terms Db{t) are changing in time and so there is a and a ti G (0, T] for which 
\Dhi{t)\ = ||Z)(t)||oo for t G [0,ti].By (4.19) we have that for this 6^ and any t G [0,ti] 

Dbi{s)Rbi ■ D{s)ds < I SbADb^{s)\\\D{s)\\^ds = [ SbiDbi{sfds, 
Jo Jo 

which, after integrating (4.21), yields 

WmWlo = Dbi{tf < D,i{Of + 2 f D,^{s)r^. • D{s)ds < p(0)||L + 2KM f \\D{s)\\1ds, 

Jo Jo 

where the final inequality makes use of (4.20). An application of Gronwall's inequality now gives us that for 

t G [0,ti] 

To complete the proof, continue this process for i > 2 by choosing the U for which is maximal 

on the time interval ti — U-i. We must have limj_^oo ti = T because (i) there are a finite number of time 
varying Z)fe(f)'s and [ii] each Dbit) is differentiable. After taking square roots we find sup^<ji ||Z)(t)||oo < 
11^(0) II ooc^*^"^ < Cie^^'^, which is equivalent to (4.9). 

We now turn our attention to showing (4.10), which we show in a similar manner. There is a C2 > such 
that for all x G and A:, £ G [1, . . . , M], 

\Dm{0,x)\ = V^\f{x + Ui/V + Uk/V) - f{x + i^e/V) - f{x + Uk/V) + f{x)\ < C2. 

Another tedious reordering of terms, which makes use of (4.18), shows that Dke{t, x) satisfies 

dtDkeit, x) = Y^ Aj{x)V[Dki{t, x + Uj/V) - Dki{t, x)] + ^{AJ{x + ue/V) - AJ {x))VDkj{t, x + ue/V) 
j j 

+ Y.{Aj{x + Uk/V) - Aj{x))VDje{t, x + Uk/V) + gM{t, x), 
j 

where 

9keit, x) = [^]{^ + ^dV + ^k/V) - Ajix + u,/V) - AJ{x + Uk/V) + AJ{x)] 

j 

X Dj{t,x + i^i/V + Uk/V). 
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By (i) the fact that the second derivative of Aj is uniformly (in j and x) bounded and (m) the bound (4.9), 
the absolute value of the last term is uniformly (in t < T, x, k, and i) bounded by some C3 > 0. 

As we did for both v(t, x) and Dk{t, x), we change perspective by viewing the above as a linear system 
with state space {k, ^, x} G [1, • • • , M] x [1, . . . , M] x = T2, where we again put an ordering on r2 
and consider MF"^ defined similarly to MF^ . Also similarly to before, we note that only a finite number of the 
Dk £{t, x) are changing in time. For b = {k, £, x} G we see that Di,{t) satisfies 

(4.22) D',,{t) = -SbVDbit) + VRb ■ D{t) + n ■ D{t) + gb{t), 

where Db{t), D{t), Sb, Rb, and rb are defined similarly as before and where we retain the necessary inequal- 
ities: for V £ MF^ 



\Rb ■ V 

(4.23) 



< ShWvl 



b2,b3+i^j/V} 

j 

In ■ v\ < 2i^M||u||oo. 

The rest of the proof is similar to the proof that the Dk{t, x) axe uniformly bounded. There is a 6^ G r2 and 
atiG (0, T] for which \Dbi (t) \ = \\D{t)\\oo for all t £ [0, h]. Taking the derivative of Dbi (t)^ while using 
(4.22), integrating, and using the bounds (4.23), we have that for this and any t £ [0, ti] 

Dbi{tf = Dbi{0f + 2 [ gbi{s)Dbi{s)ds-2 [ SbiVDbi{sfds 
Jo Jo 

+ 2 1 VRbi ■ D{s)Dbi{s)ds + 2 rbi ■ D{s)Dbi{s)ds 
Jo Jo 

<Db{0f + 2C3t + {4KM + 2C3) [ \\D{s)\\lds, 

Jo 

where we used the inequality x < l + x^ on the term D^i (s) in the first integral above. Therefore, for t < ti 

\\D{t)\\l < \\Dml + 2Cst + {4KM + 2C3) f \\D{s)\\lds. 

Jo 

We continue now by choosing a 62 G r2 suchthat |Z);,2(^)I = ||-D(t)||oo for allt G [ti, t2), with ii <t2<T. 
By similar arguments as above we have that for t G [^1,^2] 

\mt)f < \\Dih)\\l + 2Cs{t - h) + {AKM + 2C3) f \\D{s)\\1ds 

Ju 



<||M0)||^ + 2C3t + (4™ + 2C3) / \\Ks)\\i,ds. 







Continuing in this manner shows that the above inequality holds for all t G [0,T] and so a Gronwall 
inequality gives us that for all t < T 

2 </^llnrmi|2 I \jAKM+2Ci)T 



which, after taking square roots, is equivalent to (4. 10). □ 
Proof. (Of Lemma 4.4) By (4.18), we have that for any A; G [1, • • . , M] and x el7 

Dk{t,x) = Dk{^,x) + y^A]{x) [ Dkj{s,x)ds + S2{Aj{x + iyk/V)-Aj{x))vf D,{s,x + Uk/v)ds. 

J ^ Jo 

The proof is now immediate in light of Lemma 4.3. □ 
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Fig 5.1. Relative frequency of X{1) from 200, 000 sample paths constructed using (i) Gillespie's algorithm, blue line V marker, 
(ii) Euler tau-leaping, green line, O marker, and {Hi) midpoint tau-leaping, red line, * marker. The approximated distribution 
generated via midpoint tau-leaping is clearly closer to the exact distribution than that of Euler tau-leaping. 



5. Examples. 

Example 5.1. Consider the case of an irreversible isomerization of one molecule into another. We 
denote by A the molecule undergoing the isomerization and B the target molecule. We assume that the rate 
constant associated with this reaction is 1. The pictorial representation for this system is simply 

a\b. 

Letting X{t) denote the number of A molecules at time t > 0, X{t) satisfies 

X{t) = X(0) -^(yj^ X{s)ds^ . 

Supposing that we start with V = 10,000 molecules, we approximate the distribution of X{1) using 
200, 000 sample paths constructed using the Gillespie algorithm, which produces statistically exact sample 
paths, Euler tau-leaping with a step-size of 1/20, and midpoint tau-leaping with a step-size of 1/20. Note 
that in this case 1/20 = l/V'^"^^, and so /? = .325. The computational results are presented in Figure 5.1, 
which demonstrate the stronger convergence rate of midpoint tau-leaping as compared to Euler tau-leaping. 

It is simple to show that X{1) is a binomial(n, p) random variable with parameters n = 10,000 and 
p = 1/e. Therefore, EX(1) = 10000/e 3678.8. The estimated means produced from the 200,000 sample 
paths of Euler tau-leaping and midpoint tau-leaping were 3585.4 and 3681.4, respectively. Solving for £{t) 
of (3.15) for this example yields £{t) = (l/2)e^*f. Theorem 4.1 therefore estimates that Euler tau-leaping 
should produce a mean (l/2)e-il0000i--325 ~ 92.2 smaller than the actual mean, which is in agreement 
with 3678.8 - 3585.4 = 93.4. Solving for £i{t) of (3.26) for this example yields £i{t) = (l/6)te-*. The- 
orem 4.5 therefore estimates that midpoint tau-leaping should produce a mean (l/6)e~^10000^^^*'^^^ = 
4.62 smaller than the actual mean, which is in agreement with 3678.8 — 3681.4 = —2.6. 



imsart-aap ver. 2009/02/27 file: AndGangKurtzRevision.tex date: July 9, 2010 



26 



ANDERSON, GANGULY, AND KURTZ 




1000 
Prey 



Fig 5.2. Oscillations in a predator-prey model. In the left image we see the numbers of predators versus the number of prey for 
a single realization of the system 5.1. In the right image we see the time-series of the numbers of predators and prey for a single 
realization of 5.1. 



Example 5.2. We now consider a simple Lotka-Volterra predator-prey model. Letting A and B repre- 
sent the prey and predators, respectively, in a given environment we suppose {i) prey reproduce at a certain 
rate, {ii) interactions between predators and prey benefit the predator while hurting the prey, and {in) preda- 
tors die at a certain rate. One possible model for this system is 



2A, 



B 



where a choice of rate constants has been made. Letting X{t) G Z?,q be such that Xi{t) and X2{t) represent 

the numbers of prey and predators at time t > 0, respectively, X{t) satisfies 

(5.1) 



X{t) = X{0)+Yi{ / 2Xi{s)ds 



+Y2{ / m2Xi{s)X2is)ds 



1 



+^3 / 2X2is)ds 



We take X{0) = [1000, 1000]^, and so V = 1000 for our model. Lotka-Volterra models are famous for 
producing periodic solutions; this behavior is demonstrated in Figure 5.2. 

We approximate the distribution of X2(10) using 30, 000 sample paths constructed using the Gillespie 
algorithm, Euler tau-leaping with a step-size of 1/20, and midpoint tau-leaping with a step-size of 1/20. 
Note that in this case 1/20 = l/V'^^"^, and so /3 = .434. The computational results are presented in Figure 
5.3, which again demonstrate the stronger convergence rate of midpoint tau-leaping as compared to Euler 
tau-leaping. 

APPENDIX A 

Lemma A.l. Let M be a {J^t}-tna''tingale, R be bounded and {J- 1\ -adapted, and let h > 0. Then for 

r/(t) = [t/h]h, 



M{t)= / Ror]{s){M{s)-Mor]{s))ds + Rori{t){M{t)-Mor]{t)){r]{t) + h-t) 
Jo 

is an {J^t} -martingale and 

(A.1) lM]t= [ {Ror]{r)f{ri{r) + h-rfd[M]r. 
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Fig 5.3. Relative frequency of X2{10) from 30, 000 sample paths constructed using (i) Gillespie's algorithm, blue line V marker, 
(ii) Euler tau-leaping, green line, O marker, and (Hi) midpoint tau-leaping, red line, * marker. The approximated distribution 
generated via midpoint tau-leaping is clearly closer to the exact distribution than that of Euler tau-leaping. 

If M is -valued and R is Wi"^^'^ -valued, then the quadratic covariation matrix is 

[M]t= [ {7]{r) + h-rfRori{r))d[M]rR^ori{r). 
Jo 

Proof. Fort <T -h, 

lK[M{T)\Tt] = IE[ / Ro ri{s){M{s) -Mo ri{s))ds\Tt] 
Jo 

pit)+h 

= E[ Ro7]{s){M{s)-Mo7]{,s))ds\Tt] 
Jo 

= / Rorj{s){M{s)-Mori{s))ds + Rori{t){M{t)-Mori{t)){r]{t) + h-t). 
Jo 

The case ofT — /i<t<Tis similar. [M] is just the quadratic variation of the second term on the right, 
and noting that M is continuous at t = /c/i for all /c = 0, 1, 2 . . ., (A. 1) follows. □ 

For completeness we include a statement of the martingale central limit theorem (see [6] for more details). 

Lemma A. 2. Let {M„} be a sequence ofW^-valued martingales with M„(0) = 0. Suppose 

lim E[sup|M„(s) - M„(s-)|] = 

s<t 

and 

[MlMi\t^C,j{t) 
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for all t > where C = {(cij)) is deterministic and continuous. Then M„ =^ M, where M is Gaussian 
with independent increments and'E[M{t)M{t)'^] = C{t). 
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